library("dummies")
## dummies-1.5.6 provided by Decision Patterns
library("AER")
## Loading required package: car
## Loading required package: carData
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
library("plotly")
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library('RColorBrewer')
library("data.table")
library("mlogit")
## Loading required package: Formula
library("gmnl")
## Loading required package: maxLik
## Loading required package: miscTools
##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/
library("data.table")
#import and clean data
setwd("/Users/liutianyao/Desktop/Pricing Analytics/HW2")
kiwi_bubbles = read.csv('/Users/liutianyao/Desktop/Pricing Analytics/HW2/kiwi_bubbles_P2.csv')
kiwi_bubbles=kiwi_bubbles[!(kiwi_bubbles$price.KB==99),]
kiwi_bubbles=kiwi_bubbles[!(kiwi_bubbles$price.KR==99),]
kiwi_bubbles=kiwi_bubbles[!(kiwi_bubbles$price.MB==99),]
###Q3 (2)###
mlogitdata=mlogit.data(kiwi_bubbles,id="id",varying=4:7,choice="choice",shape="wide")
#Run MLE.
mle= gmnl(choice ~ price, data = mlogitdata)
summary(mle)
##
## Model estimated on: Wed Feb 12 22:59:08 2020
##
## Call:
## gmnl(formula = choice ~ price, data = mlogitdata, method = "nr")
##
## Frequencies of categories:
##
## 0 KB KR MB
## 0.41564 0.18035 0.20039 0.20362
##
## The estimation took: 0h:0m:0s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## KB:(intercept) 4.25316 0.32821 12.959 < 2.2e-16 ***
## KR:(intercept) 4.36240 0.32945 13.241 < 2.2e-16 ***
## MB:(intercept) 4.20440 0.31331 13.419 < 2.2e-16 ***
## price -3.73793 0.23671 -15.791 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -1909
## Number of observations: 1547
## Number of iterations: 4
## Exit of MLE: gradient close to zero
#demand function
demand=function(priceKB,priceKR,para){
probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(cbind(probKB,probKR))
}
#Unit cost
uc=0.5
#Profit function
profit=function(priceKB,priceKR,para){
profitKB=demand(priceKB,priceKR,para)[,1]*(priceKB-uc)*1000
profitKR=demand(priceKB,priceKR,para)[,2]*(priceKR-uc)*1000
return(cbind(profitKB,profitKR))
}
#define para and pricespace
para=mle$coefficients
priceMB=1.43
aux=seq(1,3,0.1)
pricespace=expand.grid(aux,aux,1.43)
#profit
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
profitmat[i]=sum(profit(pricespace[i,1],pricespace[i,2],para))
}
#Draw figure
xaxis=list(title="P^{KB}")
yaxis=list(autorange = "reversed",title="P^{KR}")
zaxis=list(title="Profit")
p=plot_ly(x=pricespace[,1],y=pricespace[,2],z=as.numeric(profitmat),
type="scatter3d",mode="markers",
marker = list(color = as.numeric(profitmat), colorscale = c('#FFE1A1', '#683531'), showscale = TRUE))%>%
layout(scene=list(xaxis=xaxis,yaxis=yaxis,zaxis=zaxis))%>%
config(mathjax = 'cdn')
p
#optimal price
pricespace[profitmat==max(profitmat)]
## [1] 1.20 1.20 1.43
#so the optimal price for KB and KR are 1.2 and 1.2
###Q3 (1)###
#calculate average prices
#calculate own price elasticity
KBavg = mean(kiwi_bubbles$price.KB)
KRavg = mean(kiwi_bubbles$price.KR)
MBavg = mean(kiwi_bubbles$price.MB)
probKB=exp(para[1]+para[4]*KBavg)/(1+exp(para[1]+para[4]*KBavg)+exp(para[2]+para[4]*KRavg)+exp(para[3]+para[4]*MBavg))
probKR=exp(para[2]+para[4]*KRavg)/(1+exp(para[1]+para[4]*KBavg)+exp(para[2]+para[4]*KRavg)+exp(para[3]+para[4]*MBavg))
probMB=exp(para[3]+para[4]*MBavg)/(1+exp(para[1]+para[4]*KBavg)+exp(para[2]+para[4]*KRavg)+exp(para[3]+para[4]*MBavg))
elaKB = (-para[4])*KBavg*(1-probKB) #4.257847
elaKR = (-para[4])*KRavg*(1-probKR) #4.13127
elaMB = (-para[4])*MBavg*(1-probMB) #4.069547
#calculate cross price elasticity
#KR on KB
elaKRKB = (-para[4])*KRavg*probKR #1.019923
#KR on MB
elaKRMB = (-para[4])*KRavg*probKR #1.019923
#KB on KR
elaKBKR = (-para[4])*KBavg*probKB #0.9054743
#KB on MB
elaKBMB = (-para[4])*KBavg*probKB #0.9054743
#MB on KB
elaMBKB = (-para[4])*MBavg*probMB #0.9601564
#MB on KR
elaMBKR = (-para[4])*MBavg*probMB #0.9601564
###Q4###
#Estimate single segment logit as a point of comparison
mle_noseg= gmnl(choice ~ price, data = mlogitdata)
summary(mle_noseg)
##
## Model estimated on: Wed Feb 12 22:59:08 2020
##
## Call:
## gmnl(formula = choice ~ price, data = mlogitdata, method = "nr")
##
## Frequencies of categories:
##
## 0 KB KR MB
## 0.41564 0.18035 0.20039 0.20362
##
## The estimation took: 0h:0m:0s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## KB:(intercept) 4.25316 0.32821 12.959 < 2.2e-16 ***
## KR:(intercept) 4.36240 0.32945 13.241 < 2.2e-16 ***
## MB:(intercept) 4.20440 0.31331 13.419 < 2.2e-16 ***
## price -3.73793 0.23671 -15.791 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -1909
## Number of observations: 1547
## Number of iterations: 4
## Exit of MLE: gradient close to zero
coef_noseg=mle_noseg$coefficients
#Load demographic data
demo=fread("demo_P2.csv",stringsAsFactors = F)
#Number of individuals
N = 283
#Clustering
demo_cluster = kmeans(x=demo[, 2:18], centers = 6, nstart = 1000)
# now combine cluster identity into the raw data
cluster_id = data.frame(id = demo$id)
cluster_id$cluster = demo_cluster$cluster
data = merge(kiwi_bubbles, cluster_id, by = "id", all.x = T)
# for those who don't fit in any cluster, group them into one additional cluster
data$cluster[is.na(data$cluster)] = 7
# segment share
seg.share = c( table(demo_cluster$cluster),N - sum(table(demo_cluster$cluster))) / N
# just store the coefficients (you can store many other things)
coef.est = data.frame(segment = 1:7, intercept.KB = NA, intercept.KR = NA,
intercept.MB = NA, price.coef = NA)
#Write a for-loop.
for (seg in 1:7) {
# During each loop, pick subset of data of consumers from each segment.
data.sub = subset(data, cluster == seg)
#Using that data, the rest remains the same.
mlogitdata=mlogit.data(data.sub,id="id",varying=4:7,choice="choice",shape="wide")
#Run MLE.
mle= gmnl(choice ~ price, data = mlogitdata)
mle
#Store the outcome in the coef.est matrix.
coef.est[seg, 2:5] = mle$coefficients
}
#demand functions
demandKB=function(priceKB,priceKR,priceMB,para){
probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(probKB)
}
demandKR=function(priceKB,priceKR,priceMB,para){
probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(probKR)
}
demandMB=function(priceKB,priceKR,priceMB,para){
probMB=exp(para[3]+para[4]*priceMB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(probMB)
}
pricespace=seq(0,2.0,0.01)
plot(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[1,2:5])),type='l',xlab='Prices',
ylab='Probability of purchase',col="blue",lwd=20*seg.share[1],
cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5, ylim = c(0,1))
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[5,2:5])),col="blue",lwd=20*seg.share[5])
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[4,2:5])),col="blue",lwd=20*seg.share[4])
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[3,2:5])),col="blue",lwd=20*seg.share[3])
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[6,2:5])),col="blue",lwd=20*seg.share[6])
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[2,2:5])),col="blue",lwd=20*seg.share[2])
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[7,2:5])),col="blue",lwd=20*seg.share[7])
#calculate elasticity
#Calculate profit
#choice KB
agg_choiceKB=function(priceKB,priceKR,priceMB) {
agg_choiceKB=seg.share[1]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
seg.share[2]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
seg.share[3]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
seg.share[4]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
seg.share[5]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
seg.share[6]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
seg.share[7]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))
return(agg_choiceKB)
}
#choice KR
agg_choiceKR=function(priceKB,priceKR,priceMB) {
agg_choiceKR=seg.share[1]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
seg.share[2]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
seg.share[3]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
seg.share[4]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
seg.share[5]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
seg.share[6]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
seg.share[7]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))
return(agg_choiceKR)
}
#choice MB
agg_choiceMB=function(priceKB,priceKR,priceMB) {
agg_choiceMB=seg.share[1]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
seg.share[2]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
seg.share[3]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
seg.share[4]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
seg.share[5]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
seg.share[6]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
seg.share[7]*demandMB(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))
return(agg_choiceMB)
}
#aggregate probability of KB
plot(pricespace,agg_choiceKB(pricespace,mean(data$price.KR),mean(data$price.MB)),type='l',xlab='Prices',
ylab='Choice probability',col="blue",lwd=2
,cex=2,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
lines(pricespace,demandKB(pricespace,mean(data$price.KR),mean(data$price.MB),coef_noseg),col="orange",lwd=20*seg.share[4],lty=2)
legend(1.1, 0.9, legend=c("With segmentation", "Without segmentation"),
col=c( "blue","orange"), lty=1, cex=1.1)
#KB own elasticity -4.286751
old = agg_choiceKB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new = agg_choiceKB(mean(data$price.KB)*1.01, mean(data$price.KR),mean(data$price.MB))
segelaKB = ((new-old)/old)/0.01
#KR own elasticity -3.495968
old1= agg_choiceKR(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new1= agg_choiceKR(mean(data$price.KB), mean(data$price.KR)*1.01,mean(data$price.MB))
segelaKR = ((new1-old1)/old1)/0.01
#MB own elasticity -4.078755
old2= agg_choiceMB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new2= agg_choiceMB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB)*1.01)
segelaMB = ((new2-old2)/old2)/0.01
#KR-KB cross elasticity 1.016185
old3= agg_choiceKB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new3= agg_choiceKB(mean(data$price.KB), mean(data$price.KR)*1.01,mean(data$price.MB))
segelaKRKB = ((new3-old3)/old3)/0.01
#MB-KB cross elasticity 1.01548
old4= agg_choiceKB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new4= agg_choiceKB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB)*1.01)
segelaMBKB = ((new4-old4)/old4)/0.01
#KB-KR cross elasticity 0.7162632
old5= agg_choiceKR(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new5= agg_choiceKR(mean(data$price.KB)*1.01, mean(data$price.KR),mean(data$price.MB))
segelaKBKR = ((new5-old5)/old5)/0.01
#MB-KR cross elasticity 0.8204181
old6= agg_choiceKR(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new6= agg_choiceKR(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB)*1.01)
segelaMBKR = ((new6-old6)/old6)/0.01
#KB-MB cross elasticity 0.8127758
old7= agg_choiceMB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new7= agg_choiceMB(mean(data$price.KB)*1.01, mean(data$price.KR),mean(data$price.MB))
segelaKBMB = ((new7-old7)/old7)/0.01
#KR-MB cross elasticity 0.9320796
old8= agg_choiceMB(mean(data$price.KB), mean(data$price.KR),mean(data$price.MB))
new8= agg_choiceMB(mean(data$price.KB), mean(data$price.KR)*1.01,mean(data$price.MB))
segelaKRMB = ((new7-old7)/old7)/0.01
#the most elastic one is KR-KB, the least elastic one is KB-KR
#which means they are close substitutes, but customers are more loyal to KB than KR with a higher price
#the second smallest elasticity is KB-MB, which is good for us because raising KB's
#price will not increase MB's demand so much,this reflect that KB's customers are less
#price-sensitive than MB
##customer segmentation and their preferences
#KB-KR comparison
plot(coef.est[1,2]-coef.est[1,3],coef.est[1,5],cex=20*seg.share[1],xlim=c(-3,3),ylim=c(-9,-1.5),
col = "chocolate",pch=16,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
xlab="beta_0^KB-beta_0^KR",ylab=("beta_1"))
points(coef.est[2,2]-coef.est[2,3],coef.est[2,5],cex=20*seg.share[2],col = "yellow",pch=16)
points(coef.est[3,2]-coef.est[3,3],coef.est[3,5],cex=20*seg.share[3],col = "skyblue",pch=16)
points(coef.est[4,2]-coef.est[4,3],coef.est[4,5],cex=20*seg.share[4],col = "darkgreen",pch=16)
points(coef.est[5,2]-coef.est[5,3],coef.est[5,5],cex=20*seg.share[5],col = "pink",pch=16)
points(coef.est[6,2]-coef.est[6,3],coef.est[6,5],cex=20*seg.share[6],col = "purple",pch=16)
points(coef.est[7,2]-coef.est[7,3],coef.est[7,5],cex=20*seg.share[7],col = "black",pch=16)
#KB-MB comparison
plot(coef.est[1,2]-coef.est[1,4],coef.est[1,5],cex=20*seg.share[1],xlim=c(-3,3),ylim=c(-9,-1.5),
col = "chocolate",pch=16,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
xlab="beta_0^KB-beta_0^MB",ylab=("beta_1"))
points(coef.est[2,2]-coef.est[2,4],coef.est[2,5],cex=20*seg.share[2],col = "yellow",pch=16)
points(coef.est[3,2]-coef.est[3,4],coef.est[3,5],cex=20*seg.share[3],col = "skyblue",pch=16)
points(coef.est[4,2]-coef.est[4,4],coef.est[4,5],cex=20*seg.share[4],col = "darkgreen",pch=16)
points(coef.est[5,2]-coef.est[5,4],coef.est[5,5],cex=20*seg.share[5],col = "pink",pch=16)
points(coef.est[6,2]-coef.est[6,4],coef.est[6,5],cex=20*seg.share[6],col = "purple",pch=16)
points(coef.est[7,2]-coef.est[7,4],coef.est[7,5],cex=20*seg.share[7],col = "black",pch=16)
#KR-MB comparison
plot(coef.est[1,3]-coef.est[1,4],coef.est[1,5],cex=20*seg.share[1],xlim=c(-3,3),ylim=c(-9,-1.5),
col = "chocolate",pch=16,cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5,
xlab="beta_0^KR-beta_0^MB",ylab=("beta_1"))
points(coef.est[2,3]-coef.est[2,4],coef.est[2,5],cex=20*seg.share[2],col = "yellow",pch=16)
points(coef.est[3,3]-coef.est[3,4],coef.est[3,5],cex=20*seg.share[3],col = "skyblue",pch=16)
points(coef.est[4,3]-coef.est[4,4],coef.est[4,5],cex=20*seg.share[4],col = "darkgreen",pch=16)
points(coef.est[5,3]-coef.est[5,4],coef.est[5,5],cex=20*seg.share[5],col = "pink",pch=16)
points(coef.est[6,3]-coef.est[6,4],coef.est[6,5],cex=20*seg.share[6],col = "purple",pch=16)
points(coef.est[7,3]-coef.est[7,4],coef.est[7,5],cex=20*seg.share[7],col = "black",pch=16)
#calculate each segment's profit for KB
pricespace0 = seq(0.1, 2, 0.01)
profit1=1000*seg.share[1]*demandKB(pricespace0,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[1,2:5]))*(pricespace0-uc)
profit2=1000*seg.share[2]*demandKB(pricespace0,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[2,2:5]))*(pricespace0-uc)
profit3=1000*seg.share[3]*demandKB(pricespace0,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[3,2:5]))*(pricespace0-uc)
profit4=1000*seg.share[4]*demandKB(pricespace0,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[4,2:5]))*(pricespace0-uc)
profit5=1000*seg.share[5]*demandKB(pricespace0,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[5,2:5]))*(pricespace0-uc)
profit6=1000*seg.share[6]*demandKB(pricespace0,mean(data$price.KR),mean(data$price.MB),as.numeric(coef.est[6,2:5]))*(pricespace0-uc)
pricespace0[profit1==max(profit1)]##0.99 and profit = 56.41933
## [1] 0.94
pricespace0[profit2==max(profit2)]##0.98 and profit = 21.62776
## [1] 0.99
pricespace0[profit3==max(profit3)]##1.05 and profit = 37.82587
## [1] 1.01
pricespace0[profit4==max(profit4)]##0.99 and profit = 29.18509
## [1] 0.99
pricespace0[profit5==max(profit5)]##1.01 and profit = 44.77759
## [1] 0.98
pricespace0[profit6==max(profit6)]##0.94 and profit = 25.82988
## [1] 1.05
# So we (KB) should target segment1 and segment 5
#Do not launch KB: if we estimate new coefficients#
kiwi_bubbles1=kiwi_bubbles[!(kiwi_bubbles$choice=='KB'),-5]
mlogitdata1=mlogit.data(kiwi_bubbles1,id="id",varying=4:6,choice="choice",shape="wide")
#Run MLE.
mle1= gmnl(choice ~ price, data = mlogitdata1)
summary(mle1)
##
## Model estimated on: Wed Feb 12 22:59:10 2020
##
## Call:
## gmnl(formula = choice ~ price, data = mlogitdata1, method = "nr")
##
## Frequencies of categories:
##
## 0 KR MB
## 0.50710 0.24448 0.24842
##
## The estimation took: 0h:0m:0s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## KR:(intercept) 5.20020 0.45829 11.347 < 2.2e-16 ***
## MB:(intercept) 4.97626 0.43171 11.527 < 2.2e-16 ***
## price -4.33292 0.33009 -13.127 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -1210.5
## Number of observations: 1268
## Number of iterations: 4
## Exit of MLE: gradient close to zero
para1 = mle1$coef
data1 = merge(kiwi_bubbles1, cluster_id, by = "id", all.x = T)
# for those who don't fit in any cluster, group them into one additional cluster
data1$cluster[is.na(data1$cluster)] = 7
# segment share
seg.share1 = c( table(demo_cluster$cluster),N - sum(table(demo_cluster$cluster))) / N
# just store the coefficients (you can store many other things)
coef.est1 = data.frame(segment = 1:7, intercept.KR = NA,
intercept.MB = NA, price.coef = NA)
#construct coef
for (seg in 1:7) {
# During each loop, pick subset of data of consumers from each segment.
data.sub1 = subset(data1, cluster == seg)
#Using that data, the rest remains the same.
mlogitdata1=mlogit.data(data.sub1,id="id",varying=4:6,choice="choice",shape="wide")
#Run MLE.
mle2= gmnl(choice ~ price, data = mlogitdata1)
mle2
#Store the outcome in the coef.est matrix.
coef.est1[seg, 2:4] = mle2$coefficients
}
#demand function
demandKR1=function(priceKR,priceMB,para1){
probKR1=exp(para1[1]+para1[3]*priceKR)/(1+exp(para1[1]+para1[3]*priceKR)+exp(para1[2]+para1[3]*priceMB))
return(probKR1)
}
#aggregated choice probability
agg_choiceKR1=function(priceKR,priceMB) {
agg_choiceKR1=seg.share1[1]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[1,2:4]))+
seg.share1[2]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[2,2:4]))+
seg.share1[3]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[3,2:4]))+
seg.share1[4]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[4,2:4]))+
seg.share1[5]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[5,2:4]))+
seg.share1[6]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[6,2:4]))+
seg.share1[7]*demandKR1(priceKR,priceMB,as.numeric(coef.est1[7,2:4]))
return(agg_choiceKR1)
}
#profit=1000*agg_choiceKR1(pricespace,1.43)*(pricespace-uc)
pricing1 = data.frame(pricespace = seq(0,2,0.01), profit=NA)
pricing1$profit = 1000*agg_choiceKR1(pricing1$pricespace,1.43)*(pricing1$pricespace-uc)
pricing1[pricing1$profit == max(pricing1$profit),] #maximized profit is 317.4735 when price is 1.06
## pricespace profit
## 107 1.06 317.4735
#Do not launch KB: if we use existing coefficients#
demandKR2=function(priceKR,priceMB,para){
probKR2=exp(para[1]+para[3]*priceKR)/(1+exp(para[1]+para[3]*priceKR)+exp(para[2]+para[3]*priceMB))
return(probKR2)
}
#aggregated choice probability
agg_choiceKR2=function(priceKR,priceMB) {
agg_choiceKR2=seg.share[1]*demandKR2(priceKR,priceMB,as.numeric(coef.est[1,3:5]))+
seg.share[2]*demandKR2(priceKR,priceMB,as.numeric(coef.est[2,3:5]))+
seg.share[3]*demandKR2(priceKR,priceMB,as.numeric(coef.est[3,3:5]))+
seg.share[4]*demandKR2(priceKR,priceMB,as.numeric(coef.est[4,3:5]))+
seg.share[5]*demandKR2(priceKR,priceMB,as.numeric(coef.est[5,3:5]))+
seg.share[6]*demandKR2(priceKR,priceMB,as.numeric(coef.est[6,3:5]))+
seg.share[7]*demandKR2(priceKR,priceMB,as.numeric(coef.est[7,3:5]))
return(agg_choiceKR2)
}
profit2=1000*agg_choiceKR2(pricespace,1.43)*(pricespace-uc)
pricing2 = data.frame(pricespace = seq(0,2,0.01), profit=NA)
pricing2$profit = 1000*agg_choiceKR2(pricing2$pricespace,1.43)*(pricing2$pricespace-uc)
pricing2[pricing2$profit == max(pricing2$profit),] #maximized profit is 295.7271 when price is 1.07
## pricespace profit
## 108 1.07 295.7271
#MB's profit
demandMB2=function(priceKR,priceMB,para){
probMB2=exp(para[2]+para[3]*priceMB)/(1+exp(para[1]+para[3]*priceKR)+exp(para[2]+para[3]*priceMB))
return(probMB2)
}
agg_choiceMB2=function(priceKR,priceMB) {
agg_choiceMB2=seg.share[1]*demandMB2(priceKR,priceMB,as.numeric(coef.est[1,3:5]))+
seg.share[2]*demandMB2(priceKR,priceMB,as.numeric(coef.est[2,3:5]))+
seg.share[3]*demandMB2(priceKR,priceMB,as.numeric(coef.est[3,3:5]))+
seg.share[4]*demandMB2(priceKR,priceMB,as.numeric(coef.est[4,3:5]))+
seg.share[5]*demandMB2(priceKR,priceMB,as.numeric(coef.est[5,3:5]))+
seg.share[6]*demandMB2(priceKR,priceMB,as.numeric(coef.est[6,3:5]))+
seg.share[7]*demandMB2(priceKR,priceMB,as.numeric(coef.est[7,3:5]))
return(agg_choiceMB2)
}
profitmgo1 = 1000*agg_choiceMB2(1.07,1.43)*(1.43-uc) #109.7806
#if we launch KB#
pricespaceKB = seq(0,2,0.01)
pricespaceKR = seq(0,2,0.01)
aux=seq(0,2,0.1)
pricespaceKBKR=expand.grid(aux,aux,1.43)
pricespaceKBKR[,4:6]=NA
colnames(pricespaceKBKR) = c('pricespaceKB','pricespaceKR','priceMR','profitKB','profitKR','profit')
for(i in 1:441) {
pricespaceKBKR$profitKB[i]=1000*agg_choiceKB(pricespaceKBKR$pricespaceKB[i],pricespaceKBKR$pricespaceKR[i],1.43)*(pricespaceKBKR$pricespaceKB[i]-uc)
pricespaceKBKR$profitKR[i]=1000*agg_choiceKR(pricespaceKBKR$pricespaceKR[i],pricespaceKBKR$pricespaceKR[i],1.43)*(pricespaceKBKR$pricespaceKR[i]-uc)
pricespaceKBKR$profit[i] = pricespaceKBKR$profitKB[i]+pricespaceKBKR$profitKR[i]
}
pricespaceKBKR[pricespaceKBKR$profit == max(pricespaceKBKR$profit),]
## pricespaceKB pricespaceKR priceMR profitKB profitKR profit
## 284 1 1.3 1.43 213.5893 210.8169 424.4063
#maximum profit is 424.4063 when priceKB=1, priceKR=1.3 when priceMR=1.43
profitmgo2 = 1000*agg_choiceMB(1,1.3,1.43)*(1.43-0.5) #85.3327
###Q5###
#conduct iteration to calculate Nash equilibrium
optKB = 1
optKR = 1.3
pricespace3 = seq(0, 2, 0.01)
profitmamgo = 1000*agg_choiceMB(optKB,optKR, pricespace3)*(pricespace3-uc)
max(profitmamgo) #171.3624
## [1] 171.3624
pricespace3[profitmamgo==max(profitmamgo)] #0.95
## [1] 0.95
pricespaceKBKR1=expand.grid(aux,aux,0.95)
pricespaceKBKR1[,4:6]=NA
colnames(pricespaceKBKR1) = c('pricespaceKB','pricespaceKR','priceMR','profitKB','profitKR','profit')
for(i in 1:441) {
pricespaceKBKR1$profitKB[i]=1000*agg_choiceKB(pricespaceKBKR1$pricespaceKB[i],pricespaceKBKR1$pricespaceKR[i],0.95)*(pricespaceKBKR1$pricespaceKB[i]-uc)
pricespaceKBKR1$profitKR[i]=1000*agg_choiceKR(pricespaceKBKR1$pricespaceKR[i],pricespaceKBKR1$pricespaceKR[i],0.95)*(pricespaceKBKR1$pricespaceKR[i]-uc)
pricespaceKBKR1$profit[i] = pricespaceKBKR1$profitKB[i]+pricespaceKBKR1$profitKR[i]
}
pricespaceKBKR1[pricespaceKBKR1$profit == max(pricespaceKBKR1$profit),] #priceKB = 1 and priceKR=1.2, profitKBKR = 282.6778
## pricespaceKB pricespaceKR priceMR profitKB profitKR profit
## 262 0.9 1.2 0.95 138.0391 149.6025 287.6417
#second interation
profitmamgo2 = 1000*agg_choiceMB(1,1.2, pricespace3)*(pricespace3-uc)
max(profitmamgo2) #163.5036
## [1] 163.5036
pricespace3[profitmamgo2==max(profitmamgo2)] #0.94
## [1] 0.94
pricespaceKBKR2=expand.grid(aux,aux,0.94)
pricespaceKBKR2[,4:6]=NA
colnames(pricespaceKBKR2) = c('pricespaceKB','pricespaceKR','priceMR','profitKB','profitKR','profit')
for(i in 1:441) {
pricespaceKBKR2$profitKB[i]=1000*agg_choiceKB(pricespaceKBKR2$pricespaceKB[i],pricespaceKBKR2$pricespaceKR[i],0.94)*(pricespaceKBKR2$pricespaceKB[i]-uc)
pricespaceKBKR2$profitKR[i]=1000*agg_choiceKR(pricespaceKBKR2$pricespaceKR[i],pricespaceKBKR2$pricespaceKR[i],0.94)*(pricespaceKBKR2$pricespaceKR[i]-uc)
pricespaceKBKR2$profit[i] = pricespaceKBKR2$profitKB[i]+pricespaceKBKR2$profitKR[i]
}
pricespaceKBKR2[pricespaceKBKR2$profit == max(pricespaceKBKR2$profit),] #priceKB = 1 and priceKR=1.2, profitKBKR = 278.8207
## pricespaceKB pricespaceKR priceMR profitKB profitKR profit
## 262 0.9 1.2 0.94 136.2111 147.7094 283.9205
# So, the Nash equilibrium is priceKB = 0.9, priceKR=1.2, and priceMB=0.94. Here, profitKBKR = 283.9205
#third interation
profitmamgo3 = 1000*agg_choiceMB(0.9,1.2, pricespace3)*(pricespace3-uc)
max(profitmamgo3) #144.9402
## [1] 144.9402
pricespace3[profitmamgo3==max(profitmamgo3)] #0.92
## [1] 0.92
pricespaceKBKR3=expand.grid(aux,aux,0.92)
pricespaceKBKR3[,4:6]=NA
colnames(pricespaceKBKR3) = c('pricespaceKB','pricespaceKR','priceMR','profitKB','profitKR','profit')
for(i in 1:441) {
pricespaceKBKR3$profitKB[i]=1000*agg_choiceKB(pricespaceKBKR3$pricespaceKB[i],pricespaceKBKR3$pricespaceKR[i],0.92)*(pricespaceKBKR3$pricespaceKB[i]-uc)
pricespaceKBKR3$profitKR[i]=1000*agg_choiceKR(pricespaceKBKR3$pricespaceKR[i],pricespaceKBKR3$pricespaceKR[i],0.92)*(pricespaceKBKR3$pricespaceKR[i]-uc)
pricespaceKBKR3$profit[i] = pricespaceKBKR3$profitKB[i]+pricespaceKBKR2$profitKR[i]
}
pricespaceKBKR3[pricespaceKBKR3$profit == max(pricespaceKBKR3$profit),] #priceKB = 1 and priceKR=1.2, profitKBKR = 280.217
## pricespaceKB pricespaceKR priceMR profitKB profitKR profit
## 262 0.9 1.2 0.92 132.5076 143.9094 280.217
# So, the Nash equilibrium is priceKB = 0.9, priceKR=1.2, and priceMB=0.92. Here, profitKBKR = 280.217
##without KB launch, if we work the price war, price will drop heavily and we will get less profit
# round 0
profitMB3 =1000*agg_choiceMB2(1.07,pricespace)*(pricespace-uc)
max(profitMB3)##200.1767 with price 0.99
## [1] 200.1767
pricespace[profitMB3 == max(profitMB3)]
## [1] 0.99
profitKR4 = 1000*agg_choiceKR2(pricespace,0.99)*(pricespace-uc)
max(profitKR4)##206.6413 with price 1
## [1] 206.6413
pricespace[profitKR4 == max(profitKR4)]
## [1] 1
profitMB5 = 1000*agg_choiceMB2(1,pricespace)*(pricespace-uc)
max(profitMB5)##183.5248 with price 0.97
## [1] 183.5248
pricespace[profitMB5 == max(profitMB5)]
## [1] 0.97
profitKR6 = 1000*agg_choiceKR2(pricespace,0.97)*(pricespace-uc)
max(profitKR6)##201.4765 with price 1
## [1] 201.4765
pricespace[profitKR6 == max(profitKR6)]
## [1] 1